Understanding fibrosis pathogenesis via modeling macrophage-fibroblast interplay in immune-metabolic context

Fibrosis is a progressive biological condition, leading to organ dysfunction in various clinical settings. Although fibroblasts and macrophages are known as key cellular players for fibrosis development, a comprehensive functional model that considers their interaction in the metabolic/immunologic context of fibrotic tissue has not been set up. Here we show, by transcriptome-based mathematical modeling in an in vitro system that represents macrophage-fibroblast interplay and reflects the functional effects of inflammation, hypoxia and the adaptive immune context, that irreversible fibrosis development is associated with specific combinations of metabolic and inflammatory cues. The in vitro signatures are in good alignment with transcriptomic profiles generated on laser captured glomeruli and cortical tubule-interstitial area, isolated from human transplanted kidneys with advanced stages of glomerulosclerosis and interstitial fibrosis/tubular atrophy, two clinically relevant conditions associated with organ failure in renal allografts. The model we describe here is validated on tissue based quantitative immune-phenotyping of biopsies from transplanted kidneys, demonstrating its feasibility. We conclude that the combination of in vitro and in silico modeling represents a powerful systems medicine approach to dissect fibrosis pathogenesis, applicable to specific pathological conditions, and develop coordinated targeted approaches.

. Effect of IL-4 on Mφ and Fb exposed to different metabolic and culture conditions. a-d, g-j) Volcano plots of DEGs induced by IL-4 in Mφ (a-d) and fb (g-j) when compared to resting counterparts in normoxic and single culture conditions (a, g), in hypoxic and single cultures conditions (b, h), in normoxic and cocultures conditions (c, i), in hypoxic and cocultures conditions (d, j). Number of DEGs (FDR < 0.05, threshold line: 1.3) upregulated (logFC ≥ 1) and downregulated (logFC ≤ -1) are reported in red. e, f, k) Ingenuity pathway analysis identified functional pathways with a significant positive (zscore ≥ 2, in orange) or negative (z-score ≤ -2, in violet) enrichment associated to DEGs in Mφ (e, f) and Fb (k), under same experimental conditions reported for the corresponding volcano plot.
Enriched pathways were used to generate heatmaps reported in Fig. 3e Figure S4. Transcript levels of hypoxia-responsive genes.
Hypoxia-responsive genes were investigated at different time points in single cultivated Mφ      Figure S6. Effect of cellular interplay on Mφ and Fb exposed to different metabolic and immune conditions.

a-d, h-k) Volcano plots of DEGs induced by coculture in Mφ (a-d) and
Fb (h-k) when compared to single cultivated counterparts in normoxia without immune challenge (a, h), in hypoxia without immune challenge (b, i), in normoxia with LPS+IFNγ stimulation (c, j), or in hypoxia with LPS+IFNγ (d, k).
e-g, l-m) Ingenuity pathway analysis identified functional pathways with a significant positive (zscore ≥ 2, in orange) or negative (z-score ≤ -2, in violet) enrichment associated to DEGs in Mφ (eg) and Fb (l-m), under same experimental conditions reported for the corresponding volcano plot.
Enriched pathways were used to generate heatmaps reported in Fig. 5e and 5f for Mφ and Fb, respectively.        The graphs depict the glomerular filtration rate (GFR) and the time points when consecutive biopsies were performed for the 9 patients depicted in Fig. 10c, selected to represent a large range of different underlying pathologies and patient characteristics (see Table S5). These cases  Fig. 10c show the cellular infiltrates at these time points in detail.

Multilevel analysis of transcriptomic data generated in the in vitro model.
The relative contribution of each variable and multiple combined effects was evaluated by supervised analysis at three levels of increasing complexity. In the 1 st level analysis, one by one variable was taken into account: differential gene expression analysis was assessed on normalized counts and performed in paired. For each single comparison, we defined an ad hoc design matrix taking into account replicates, phenotype, culture condition, O 2 status, and polarization information.
For each gene statistical significance of differential expression was tested using the QL F-test, only genes with a False Discovery Rate (FDR, p-value adjusted considering Benjamini-Hochberg correction) ≤ 0.05 were selected and referred to as "differentially expressed genes" (DEG). DEGs positively and negatively regulated were then shown in volcano plots. DEGs associated to each comparison and their logFC values (|logFC| ≥ 1) were used for the identification of enriched pathways using the Ingenuity Pathway Analysis software (IPA v01-13; Qiagen). Only pathways with |z-score| ≥ 2 and log 10 (p-value) ≥ 1.3 were selected. 1 st level comparisons were based on each of the following variables:

Pre-ranked Gene Set Enrichment Analysis (GSEA).
GSEA was performed in order to evaluate which in vitro signatures were significantly enriched in ex vivo data. In the first step we defined in vitro signatures: for each previous performed differential expression analysis of in vitro data, two signatures were created: an "UP regulation" signature (gene subset parameters: FDR ≤ 0.05 and logFC ≥ 1) and a "DOWN regulation" signature (FDR ≤ 0.05 and logFC ≤ -1). Only signatures with more than 15 genes were taken into account. If the length of the signature was greater than 250 genes, only first 250 (ranked by logFC) genes were considered.
Then we generated four ranked lists (defined on logFC values) referred to the general comparison of pathological and control samples and to the same comparison but specific for the three distinct anatomical regions (interstitium, glomeruli surrounding area and glomeruli). Finally, pre-ranked GSEA was performed for each ranked list applying 1000 permutations, weighted enrichment statistic method; seed set for permutation was 149.
The overlap analysis was performed with Gene Overlap Bioconductor package, which evaluates pvalue, odds-ratio and Jaccard index of in vitro and ex vivo signatures overlaps.

Percentage explained analysis.
Enriched in vitro signatures and overlap significance for each ex vivo region are known.
Considering one by one region, only in vitro signatures respecting specific filters were taken into account (for GSEA: Normalized Enrichment score greater than 1, nominal p-value lower than 0.01; for Gene Overlap: p-value < 0.05 and odds-ratio > 1).
Let signature A, B, C etc. be in vitro signatures (with AA as last signature), let signature α be the ex vivo signature of the selected region, let signature M be the Metascape signature, K the number of the enriched Metascape signatures and let n assume the meaning of "number of", considering signatures as defined in "Gene overlap" section, we determined: a) percentage of ex vivo genes explained by in vitro model (with genes defined as union of all genes of the selected in vitro signatures). The results are represented in the pie chart. Specific pathway analysis.
Considering glomeruli, interstitium, and glomeruli surrounding area, genes not shared were marked as region-specific genes. Signatures obtained as just described were been input data in new pathway analyses. Pathway resulted were subsets of previous pathway analyses (see the Methods section) and were highlighted with asterisks.

Real-time quantitative PCR (qPCR).
Total mRNA was extracted from human macrophages and fibroblasts using DirectZOL RNA

Robustness analysis of model's dynamic behaviour.
To assess the robustness of the model's qualitative behaviours, we our analysis with respect to variations all free parameter in the model.

Bifurcation analysis.
As a first step, we performed a bifurcation study by fixing all parameter values and varying the values of each parameter individually (Fig. S8). We found that all parameters are bifurcating, except for the parameter α, which only affects the value of the fibrotic steady state. With respect to all other parameters, we found that two bifurcations take place: for a considerable parameter range, only the healthy steady state exists and is globally stable. At certain parameter values, a saddle-node bifurcation occurs which creates a stable fibrotic steady state, and a saddle point. In this parameter range, the system is bistable and the initial conditions dictate whether the system ends up in the fibrotic or healthy steady states. Further change of the parameter values, a second bifurcation happens when the saddle node and healthy steady states collide, a transcritical bifurcation. Thus, the healthy steady state loses stability, and the fibrotic steady state becomes the only biologically relevant stable steady state. No other bifurcations occur. Since the system is two-dimensional, no chaotic behaviour or hidden attractors are possible. Thus, the only possible long-time behaviour according to the model is observing the fibrotic or healthy steady states. Furthermore, we confirmed this long-time dynamic behaviour through numerical calculation of the ODE system solutions (Fig.   S9). For all parameters, we find that after a long time, solutions either reach the healthy or the fibrotic steady state, according to the stability criterion of the healthy steady state reported in the main text, when the initial fibroblast and macrophage concentrations are low. When the initial macrophage and fibroblast concentrations are high, we observe the bistable regime (after the saddlenode bifurcation has happened) as a small region in the parameter space where the solutions reach the fibrotic steady state even below the instability threshold of the healthy steady state. The sole noteworthy behaviour we find from simulations is visible in the H vs. y parameter space. While for all other parameters, the bistability region remains small, in this case the bistable region grows with decreasing inflammation y, which may seem counterintuitive at first. However, we should recall that the fibrotic steady state in the bistable regime can only be reached for high initial macrophage and fibroblast concentrations. If inflammation is low (or absent) then only pro-fibrotic macrophages and fibroblasts are present, which cannot decrease significantly on their own, thus leading to fibrosis. In other words, high levels of active fibroblasts in the absence of damage can only lead to fibrosis.